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SUMMARY 

Models of physical processes are often too complex to be handled directly. Approxi- 
mations of these models that are of reduced complexity but that maintain the essential 
model characteristics are often used for analysis rather than the original model. One 
way to generate these approximations is the Routh approximation technique. In this re- 
port, this technique is programmed and applied in the frequency domain to a 16th-order, 
state -variable model of the F100 engine and to a 43rd-order, transfer-function model 
of a launch vehicle boost-pump pressure regulator. The accuracies of reduced-order 
approximants obtained for both models are compared with the original model. The re- 
sults indicate that the Routh method is an excellent technique for linear system reduction 
in the frequency domain. Therefore, the frequency-domain formulation was extended to 
the time domain in order to handle the state- variable model directly. The time-domain 
formulation is derived herein and a new characterization that specifies all possible Routh 
similarity transformations is given. This characterization is computed by using 
eigenvector- eigenvalue techniques and is therefore quite accurate. The application of the 
time-domain Routh technique to the state- variable engine model is described, and 
reduced-order approximation results are given. Additionally, an optimization procedure 
that can improve the approximation accuracy by taking advantage of the transformation 
characterization is applied to the engine example. Also, the limitations of applying the 
Routh technique to the fixed-dimension, multiple- input - multiple- output reduction prob- 
lems are discussed. 


INTRODUCTION 

In control system design the system to be controlled is often represented by a com- 
plex mathematical model. Practically, this complex model may be difficult to use for 
design purposes. Additionally, the resultant control design may be too complex to im- 
plement. To eliminate these problems, model reduction methods are often used before 
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designing the control system in order to reduce the complexity of the original model 
while maintaining its important characteristics. Alternatively, model reduction methods 
could be applied to the reduction of complex control designs. In either case the model 
reduction process would be the same. 

One philosophy adopted in reduction methods in dominant-mode approximation 
(ref. 1). In this report, a dominant- mode- approximation technique, the Routh approxi- 
mation technique (refs. 2 and 3), is applied in both the frequency and time domains to 
simplify a state -variable turbofan engine model, and in the frequency domain to simplify 
a transfer-function model of a launch vehicle boost-pump pressure regulator. Hutton 
and Rabins (ref. 4) have shown the frequency-domain formulation of the Routh approxi- 
mation technique to be an excellent reduction method for many mechanical systems. 

This formulation is applicable to turbofan engine and boost-pump regulator models as 
well and is therefore a strong motivation for the work of this report. 

The objectives considered included the application of the Routh approximation meth- 
od in the frequency domain to the F100 turbofan engine. This application demonstrated 
that adequate, reduced- complexity models can be determined by the Routh technique. 
Next, the applicability of the frequency-domain approach to high-order systems was 
studied by using the boost-pump pressure regulator model. Finally, the frequency- 
domain Routh formulation was extended to the time domain. This extension allows the 
usual engine time-domain formulation to be handled directly. The time-domain Routh 
approximation formulation also allows additional flexibility in the choice of the approxi- 
mate system. By using the additional flexibility of the time-domain formulation, opti- 
mized approximations can be readily generated. The time-domain formulation was ap- 
plied to the engine model, and some results are given. Optimization results and the ap- 
plication of the Routh method to deriving fixed-dimension realizations of the F100 engine 
are also discussed. 


ROUTH APPROXIMATION 

Many physical processes can be represented as multiple- input - multiple- output 
(MIMO) dynamic systems by using linear, constant- coefficient, differential equations to 
model small variations at various operating- point conditions. A transfer-function repre- 
sentation of these equations can be used as input information in a classical or 
multivariable-frequency-response control system design. The complexity of these 
transfer functions can cloud the design process with unnecessary detail or difficult com- 
putational problems. Additionally, control systems may be made more complex by un- 
necessary model complexity. Thus, there is a strong motivation to reduce the complex- 
ity of system and control models. 
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Frequency Domain 


The Routh approximation method (refs. 2 and 3) is a dominant- mode reduction tech- 
nique. In the frequency domain the method incorporates parameters obtained from a 
Routh stability table analysis in the reduction process. These parameters, called alpha 
and beta parameters, are obtained from the coefficients of the MIMO transfer function 


W(s) = 


BjS 11-1 + . . . + B n 


a Q s n + a^s 


n-1 


+ • • . + 3 


n 


(i) 


where the a^ are scalars and the B. are l x m matrices. (All symbols are defined 
in appendix A. ) Algorithmically, the n alpha parameters are computed from the de- 
nominator of equation (1) as (ref. 4) 
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and 
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where 


aj = 0 for 


i + j > n 


The coefficients aj can be arranged as entries in a Routh table (ref. 2). This tabular 
construction suggests the name Routh approximation. 

The beta parameters are computed from both the numerator and the denominator of 
the given transfer function. Consider the single-input - single-output (SISO) transfer 
function from input q to output p implicitly given in equation (1), Wpq(s). The beta 
parameters for this SISO transfer function are 
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3 


n-j 


b rVi-i 


3 = 0,2,4, . . . , <n 


(4) 


3 


I 



and 


u n 

^7 

a 0 


"N 


K i+2 i* 
b i b j+2 


^ a j+ 2 J 


1 = 1 , 2 , 

] = 0,2,4, , . 


. . , n 
< n - i 


(5) 


where 

bj = 0 for i + j > n 

ii. 

and b f is the (qp) element of the matrix B r , where r = 1 to n. In this way, beta 
parameters for the lx m possible SISO transfer functions of equation (1) can be found. 

ii 

Again, considering the SISO transfer function Wpq(s), the computation of the k in approx- 
imant from the first k alpha and beta parameters is given recursively as 
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( 8 ) 


where 



for i + j > k 


Note that, for different output- input pairs, different approximant orders k can be 

ii ii 

chosen. The k -order approximant for the (pq) n output-input pair is written as 
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(9) 



j _,k~ 1 j k~ 2 j 

d^s + d^s + . . . + d^ 
k k-1 

C Q S + CjS + . . . + Cj^ 


Pr operties . - The Routh approximation method exhibits several very useful proper- 
ties. Each is important in the reduction problem and is described briefly here. 

Stability: If the original transfer function is asymptotically stable, the Routh approxi- 
mant, of any order, will be asymptotically stable. 

Pole-zero locations: The poles and zeros of the approximants approach the poles 
and zeros of the original function as the order of the approximation is increased. 

Impulse response energy: If 4(t) is the impulse response of W (s), and if the 
impulse response energy is defined as 



( 10 ) 


(ID 


Also, if /^(t) is the impulse response of the k -order approximant of Wjs), then 


E k+1 ‘ E k + 


^k+1 


2 a 


( 12 ) 


k+1 


Since W (s) is assumed to be asymptotically stable, a. > 0 for i=l, . . .,n and 

pcj 1 

0 < Ej < E 2 < Eg < . . . < E n = E (13) 

The ratio E v /E gives an indication of the approximation accuracy in terms of the per- 

K 

centage of total energy accounted for by the k tn approximation. 

Derivatives. - The k tn Routh approximant satisfies the following derivative con- 
dition: 

i = 0, 1, . . ., k - 1 (14) 

s=0 


— [W(s)] 


ds 


s=0 


- 0V s >] 

ds 1 


Initial and final values . - The initial and final values of the step response trajecto- 
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ries for the original system and its k -order Routh approximant are the same. This is 
true even though the approximant is a strictly proper transfer function for all p and q. 


Time Domain 


The Routh approximation method was extended to the time domain in order that Routh 
approximants could be directly determined from a state-space model formulation. The 
alternative would be to determine the pxq transfer functions for the MIMO system and 
then reduce each function individually. For a typical engine model with p = 16 and 
q = 5, this alternative would involve reducing 80 transfer functions at each operating 
point. Computationally, this is undesirable. Hutton (ref. 3) outlines a method for di- 
rectly determining the alpha and beta parameters from the state- space formulation. 
However, computationally, this method failed for the 16th- order engine example of this 
report. Thus, an alternative computational procedure was devised. It incorporates 
eigenvalue- eigenvector solutions that are well known and for which excellent computation- 
al solutions exist (ref. 5). 

Given the state- space representation S-,, 



(15) 


where x c R", u C R m , and y C R‘, define a reciprocal system, Eg, such that 

x = Ax + Bu' 
y = Cx 


a 


(16) 


where 


*■ _ i 

A = A 1 
B = -AB 
C = C 


(17) 


The system Ej is assumed to be stable and have distinct eigenvalues. For the jet en- 
gine system studied herein, this is a realistic assumption. Note that the reciprocal 
transformation of Eg is S^. This transformation will always exist since A is as- 
sumed to be stable. The transformation is required to preserve the dominant or low- 
frequency information of Sj during the reduction process (ref. 3). Now consider a 
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nonunique similarity transformation of S 2 , called the Routh transformation, such that 

at r = T r R (18) 


where R is the Routh stability matrix. That is, 

R = A R T 


(19) 


-1 -1 
10-1 

10-1 N o v 


r = 


(20) 


0 -1 

1 o 


and 


A r = diagfcr] i = 1, . • • , n 


Also, consider the modal transformation of S 2 , where 


■^m^m 


T D A 
mm 


( 21 ) 


( 22 ) 


and 


A = diagfXj] i = 1, . . . ,n (23) 

Here the Aj are eigenvalues of A and the i column of T m is the corresponding 
eigenvector. Note that T m will be unique with respect to a scaling convention and a 
specific order of the Aj’s. The scaling convention is represented by D m , an arbitrary, 
full rank, diagonal matrix with diagonal elements d. that may be complex. 

Since R is similar to A, 
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(24) 


RT D = TDA 
z z z z 


Again T z is a modal matrix for R and A, and it is unique with respect to a scaling 

convention represented by D . From equations (18), (22), and (24) and the definition 

z 


D 


= D D J 
m z 


(25) 


all possible Routh transformations can be characterized as 


T 


R 


= T DT 
m 


-1 

z 


(26) 


Recall that a fixed ordering of eigenvalues is given, 
be represented by a permutation matrix E^, where 

T R, p = E p T R 


Permutations of this ordering can 


(27) 


and Ep is the identity matrix with appropriate columns interchanged. Although this 
permutation does not affect the input-output transfer relationship of the original system, 
it will affect the input-output characteristics of the approximate system. Undoubtedly, 
one of the n! possible state permutations will yield a better approximation in some 
sense than another, but this problem is not considered herein. 

From equation (18) a Routh canonical system Eg can be written as 


where 


x R = Rx r + Gu 


y = Hx 


R 


(28) 


x 


= t r x r 


(29) 


The Routh approximation procedure starts by assuming that 

(x r ) 2 = 0 (30) 

where 
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and C R^, and |x^ C R 11 ^. The system Eg that incorporates the assumption 

of equation (30) becomes z!,, or 


*■ 1 " T R 11 R 11 T R 11 X 1 + T R n G l U 


y = h i t r u x i 


where 


'_ x _l_" 
_ x 2 _ 


and XjC R and the subscripts indicate conformable partitioning of the appropriate 
matrices. The desired approximation then is the reciprocal transformation of ZJ,, 
called Sj. 

Examination of equation (32) shows that the reduced- order system matrix, 

T R„ , R 11 T R ^ , is simply a similarity transformation of the truncated Routh stability 

matrix. Truncating an n -order Routh matrix R yields another Routh matrix of k Ln 
order R^. Since the original system is stable, the alpha parameters of R, a. (i = 1 
or n), are positive. Therefore, R^ is also stable. The reduction process implied 
in equation (32) exhibits the same pole property as the frequency-domain approach. 

Hutton (ref. 3) has shown that if the system given by equation (16) is a single- input - 
multiple -output system and T R is selected such that 


TrB = 


the elements of the rows h. of H become the beta parameters for the output considered 


h i ^li’ ^2i» ' • ' ’ 
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and the properties of the approximation described previously hold. It can be shown that 
if D is selected as D=diag[# i ], where Q = T^B, equation (33) is satisfied. 

COMPUTATIONAL ASPECTS OF ROUTH APPROXIMATION 

Computational aspects of both the frequency- and time-domain formulations are 
discussed in this section. 


Frequency Domain 

In the frequency domain the Routh approximation problem is solved by calculating 
the alpha and beta parameters for a given transfer function. Straightforward program- 
ming of equations (3), (5), and (6) as detailed in subroutines DROUTH and DTHCNV in 
appendix B gives these parameters (DROUTH) and the Routh approximants (DTHCNV). 
Also, the impulse response energy property of equation (11) was programmed in sub- 
routine DIMPLS. Implementing this Routh tabular procedure is straightforward, and 
computationally the procedure is an efficient and accurate method of reducing SISO 
transfer functions in the frequency domain. 


Time Domain 

Computationally, given the state-space formulation of equation (15), time-domain 
Routh approximation requires (1) the alpha coefficients and (2) the Routh transformation 
matrix. The alpha coefficients can be found readily from equation (3) if the system- 
characteristic- equation coefficients can be found. A program that incorporates the 
method of Daniel evski (described in ref. 6) to find the characteristic equation of the sys- 
tem matrix A was used in this report. The alpha coefficients were found by applying 
subroutine DROUTH. 

A method for computing the Routh transformation matrix of equation (18) has been 
proposed by Datta (ref. 5). This method was programmed and applied to the 16th-order 
jet engine model. Numerical results were unacceptable because of large computational 
errors, and so another technique was sought. Since accurate and efficient eigenvalue- 
eigenvector techniques are well known (for the case of distinct eigenvalues), the comput- 
ation of the Routh transformation matrix was reformulated as the solution of two 
eigenvalue- eigenvector proglems. These two problems are represented by equations 
(22) and (24) and are solved by using the computational methods of reference 7. Since 
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A and R are similar, they have the same eigenvalues. The Routh matrix R is also 
known from equation (19) since the alpha parameters have been found. 

The transformation matrices of equations (22) and (24) may have elements in the 
complex field. Thus, from equations (2) to (6), T R can in general be complex. How- 
ever, a real T^ is desired to yield a physically realizable approximation and to facili- 
tate computer computations. The matrix Tj^ can be constrained to be real by the prop- 
er selection of the elements d^ of D. To determine these elements, consider the eigen- 
vector solution of equation (24) by means of real-number computer operations. Given A 
and its eigenvalues, the modal matrix of A can be written as 


T 

mm 


= T V 
mm 


(35) 


where V„ is the block diagonal and T„^, the modified modal matrix of A, is a ma- 
in ° mm ’ ’ 

trix of real numbers. If the eigenvalues of A were all real, V„ would equal I. How- 
ever, if X^ and X^ + j were a complex conjugate pair, the i in block of V , ( V m \ , 
would be defined as \ z 1 


Likewise, for equation (24), 



T 

z, m 


= T Z V 


m 


Now from equations (26), (33), and (35) and the definition 


D B =V m DV m 


(36) 


(37) 


(38) 


Tp^ can be written 


T = T D T x 
X R 1 min B 1 z, m 


Now if Tpj is real, D^, a block diagonal matrix, must be real. This can be assured by 
selecting the elements of D, d^, such that 

(1) dj is real if is real 

(2) d 4 ,d i+1 are a complex conjugate pair if Xp X J+1 are a complex conjugate pair. 

The reciprocal transformation of equation (17) does not represent a significant in- 
crease in computations since the modal transformation is the same for A and A and 
the eigenvalues are simply reciprocals. 

The matrix inversion of equation (26) can be eliminated by noting that 
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(40) 


r p 

r = t r t 


where 

T = diag[(-l) i_1 ] 


From equation (19) then 


A r TR = R t a r T 


Now from equations (24), (26), and (40) 


T R = Wz TA R 


(41) 


(42) 


(43) 


A program that solves for T R , the Routh transformation in equation (43), was writ- 
ten. This program, called DRTTRC, not only finds T R> but also the transformed ma- 
trices of equation (28), G and H. The FORTRAN listing of DRTTRC is given in appen- 
dix B. Additionally, a program that finds the reduced system of equation (32), called 
DTYPEI, was written and is given in appendix B. 


APPLICATION TO ENGINE MODEL 

The frequency-domain formulation of the Routh approximation technique was applied 
to the state -variable model of the F100 engine at the intermediate-power operating point 
and to the transfer-function model of the boost-pump pressure regulator design. The 
time-domain formulation was applied only to the engine model. 


Frequency-Domain Applications 

As an initial experiment the frequency- domain formulation was applied to two 16th- 
order SISO transfer functions. These transfer functions represent dynamics of the tur- 
bofan engine from the input w^ (fuel flow) to the outputs N c (compressor speed) and 
Tj. (turbine inlet temperature). The coefficients for these transfer functions were calcu- 
lated and are given in table I. Table II shows the impulse response ratios for Routh ap- 
proximations of increasing order for the two transfer functions. These ratios were cal- 
culated from the alpha and beta coefficients as outlined in equation (12). From these 
ratios an acceptable approximation order can be estimated. 

The approximant order was selected by first choosing a minimum acceptable level 
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of accuracy, as defined by the impulse response energy ratio. The level chosen corres- 
ponds to a ratio of 0.81. This level of accuracy was assumed to be adequate for the 
purpose of this report. Next the order that most nearly corresponded to the selected 
ratio was chosen as the order of the Routh approximant. For the energy ratios given 
in table n, the Routh approximant orders, as found by these criteria, are 

k l ■ 3 <?1 = N c> 

k 2 = 9 (y 2 = T t ) 

The Routh approximants for N c and T^. were calculated by using the approximant 
orders and kg and the algorithm given in equations (6) to (9). The results are 
summarized in table m. Comparisons of pole locations and step responses for exact 
and approximate transfer functions were made to evaluate the adequacy of the approxi- 
mations. Table IV gives the pole locations, and figure 1 the step response comparisons. 
The pole comparison shows a good correlation between actual and approximate system 
poles. Likewise, the step responses of figure 1 indicate an excellent agreement between 
actual and approximate model representations. Based on these results it was concluded 
that the Routh approximation technique is a viable approach to reducing the complexity of 
frequency-domain models of jet engine dynamics. Also, the accuracy level selected in 
this initial example may be too stringent, based on the step response comparison. 

For comparison, the Routh approximation technique was applied to all 80 (l = 16, 
m = 5) possible transfer functions of the given engine example. The impulse response 
energy ciriterion for model order was used to determine the order for each of the ap- 
proximants. These reduced orders are given in table V. Examination of this table gives 
dynamically fast (x^, x^, e.g. ) and slow (x^, x^, x^q, e.g. ) states as well as reduced 
order. The large difference in reduced order for transfer functions with the same in- 
put indicates that modeling of the fast states as constants with respect to the slow states 
would be acceptable. 

The frequency response technique was also applied to a high-order (4 3rd order) 
model of the boost-pump pressure regulator, which represents the linearized line dy- 
namics of a liquid-oxygen supply system in a rocket engine (ref. 8). The model input 
represents a pressure error signal; the output represents an actuator piston- commanded 
velocity. The transfer function (time scaled by a factor of 100) for this system is given 
in table VI. This application tests the accuracy and the computational feasibility of the 
programs DROUTH, DTHCNV, and DIMPLS for high-order problems. Again, straight- 
forward application of the programs was all that was required to obtain Routh approxi- 
mants to the boost-pump pressure regulation system. Additionally, it was determined 
that a 43rd-order implementation of this pressure control would be impractical and that 
a study of lower order implementations would be desirable. Thus, this application is 
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of more than numerical significance. 

For the 81-percent impulse response energy criterion, a 26th-order Routh approxi- 
mant was selected. The poles of the original and approximant system are compared in 
table VIE; the zeros are compared in table VIII. Both show good comparability between 
original and approximant values. For reference the 26th- order approximant is given in 
table IX. Figure 2 compares the frequency responses for the original and approximant 
systems. The comparison is exact for both magnitude and phase for frequencies less 
than 20 hertz. The approximation is still good for higher frequencies, and the loss of 
frequency information due to the approximation is clearly shown. 


Time-Domain Application 

The time-domain formulation of the Routh application technique was applied to the 
16th-order, state- variable F100 engine model. This model represents the dynamic 
response of the F100 engine for small perturbations about the intermediate- power oper- 
ating point. The normalized system matrices are given in table X and the state, input, 
and output vector definitions in table XI. Calculating the Routh transformation matrix 
Tj^ poses the only significant computational problem in applying the time-domain for- 
mulation. Thus, as a test, the alpha parameters calculated from the characteristic 
equation of A (Routh table approach) and from the similarity transformation T^ATj^ 
are compared in table XII. Clearly, calculating Tp with DRTTRC is accomplished 
with a high degree of accuracy. Therefore, the tools exist to accurately apply the time- 
domain formulation to high-order systems. 

Examining table V shows that a fifth-order approximation for each input would 
represent an adequate tradeoff between accuracy and complexity. Reduced approximants 
were calculated by using DRTTRC and DTYPEI for this order. Step response trajec- 
tories comparing the original and approximate systems were generated. The T^. result 
indicates the degradation of accuracy obtained with a lower order model. Typical step 
responses for one input (fuel flow) and two outputs (compressor speed and turbine inlet 
temperature) are given in figure 3. 

With this time-domain formulation, fifth-order models were calculated for each 
input that exactly match those that would have been determined with a frequency-domain 
analysis. However, the approximants were determined directly from the state- variable 
model without first calculating the frequency transfer function for each input-output pair. 
All the properties of the frequency-domain approach hold for the time-domain approach 
when considering the state- variable formulation one input at a time. 
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Fixed-Dimension Approximant 


For a five- input system approximated by fifth-order models, a total reduced system 
realization may, require a state vector dimension of 25. In fact, the realization dimen- 
sion will probably vary for different numbers of inputs and orders of approximation. 
Often, however, the allowable dimension of the total reduced approximation is fixed by 
some engineering or economic constraint. Constructing a fixed-dimension realization 
of the total system from the Routh approximation can be a very difficult computational 
problem, especially if the input and output dimensions are large. Thus, the technique 
as outlined does not directly handle the fixed-dimension problem. However, if the same 
Tj^ matrix is used in the time-domain reduction process for each input, the total ap- 
proximation can be realized by a system of k** 1 order, where k is the number of alpha 
parameters retained, and the fixed-dimension approximation problem is solved. Select- 
ing the T R matrix that will give the best k -order approximation is therefore an in- 
teresting problem. 

In this regard the somewhat arbitrary selection of D can be used to good advantage. 
Many different transformations can be found very quickly once the original two eigenvec- 
tor problems have been solved. Indeed, the selection process can be automated by opti- 
mizing some function of error between the original and total approximant systems over 
the n-parameter space of D. Such an optimization scheme was tried on the example 
given in this report. Two different error functions were used for comparison. The first 
was a weighted sum of the differences in system and approximant step response energies. 
The second was a weighted sum of squares of the differences in system and approximant 
steady- state values. Significant minimization of each error function was easily achieved 
in relatively small amounts of computer time by using a conjugate direction optimization 
scheme. Thus, the general optimization procedure would appear to be a good way to 
improve the accuracy of approximants for certain systems while maintaining a fixed or- 
der of realization. However, the time-domain Routh approximation procedure, when 
constrained to yield a fixed-dimension realization in the multiple- input case, does not 
exhibit the final- value property of the single-input case. For the engine example posed, 
the significant improvement in the approximation gained by function optimization was 
overshadowed by these final-value errors for the multiple- input case. 

Based on this observation the time-domain formulation was modified to ensure that 
the final-value property would be met by a fixed-dimension, multiple- input Routh ap- 
proximant. In the original formulation the difficulty with final values in the multiple- 
input, fixed-dimension problem can be traced to the original assumption of the reduction 
process (eq. (30)). If the initial assumption were changed to 
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the approximation would force the final-value property. Unfortunately, the assumption 
of equation (44) when applied to the original system may yield an unstable approximant 
for a given stable system. This was the case for the engine example and, thus, the 
modification of the Routh procedure was rejected. 


CONCLUSIONS 

Models of physical processes are often too complex to handle directly. Approxima- 
tions of these models that are of reduced complexity but that maintain essential model 
characteristics are often used for analysis rather than the original model. One way to 
generate these approximations is the Routh approximation technique. The frequency- 
domain formulation of the Routh technique was applied to transfer-function models of an 
F100 engine and a launch vehicle boost-pump pressure regulator. Also, the Routh ap- 
proximation process was reformulated in the time domain. A new characterization of 
the nonunique Routh similarity transformation was derived that describes all possible 
Routh transformations. This characterization casts the computation of the Routh trans- 
formation into two eigenvector- eigenvalue problems that are easily solved. The applica- 
tion of the time-domain formulation to a 16th-order state-variable description of a tur- 
bofan engine was described, and the results were given. These results indicate that the 
time-domain Routh approximation technique can be valuable in reducing engine model 
complexity when dealing with the model on a single- input basis. An optimization proce- 
dure was discussed that can significantly improve the approximation in a computationally 
efficient manner by talcing advantage of the new time-domain characterization. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, May 11, 1977, 

505-05. 
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APPENDIX A 


SYMBOLS 


A reciprocal system matrix 

a 

A original system matrix 

a coefficient of characteristic equation 

B reciprocal control matrix 

B original control matrix 

B . matrix of numerator coefficients 

b element of B i 

C reciprocal output matrix 

C original output matrix 

c coefficient of characteristic equation for reduced system 

D diagonal parameter matrix 

d coefficient of numerator polynomial for reduced system 

d/ds 1 i -order differentiation operator 
E impulse response energy 

Ep permutation matrix 

G Routh control matrix 

H Routh output matrix 

A impulse response 

In number of rows 

K reduced-system order 

k Routh approximant order 

l number of outputs 

m number of inputs 

N c compressor speed 

n original- system order 

R Routh- system matrix 

s Laplace transform variable 



T similarity transformation matrix 

turbine inlet temperature 
t time, sec 

u control vector 

V block diagonal modification matrix 

W(s) matrix transfer function 

Wj engine fuel flow 

x reciprocal state vector 

£ original state vector 

y output vector 

a alpha parameter 

0 beta parameters 

r gamma matrix 

X diagonal parameter matrix 

S state- space system 

Subscripts: 

k k -order approximant 

m modal 

mm modified modal 

p output number 

q input number 

R Routh transformation 

z Routh transformed system 
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APPENDIX B 


SUBROUTINES USED IN ROUTH APPROXIMATION 


c subroutine ROUTH 

c 

C PURPOSE 

C TO CALCULATE THE ALPHA-BETA EXPANSION 

C COEFFICIENTS OF A LINEAR, YIME-INVARJANT 

C TRANSFER F LNCT ION < SI S 0 » 

C 

C HE 1HOD 

C ROUTH TABLE 

C 

SU TROUT IN E DROLTH C ?1 , 7 2 , ALPH A , R E TA , RT A , RT B , N , NH A X I 
IMPLICIT DOUBLE PRECISION tA-H,0-ZI 

D I PENSION HU 1,22(1) , ALPHA! 1 I ,BETA tl),RT ACNMAX ,1 I ,RTP(NMAX,1) 
KK=0 

DC 5 1=1, N 
DO 5 J = I , N 
RTACI ,JI=U. 

PTt (I,J1=0. 

I CONTINUE 

INITIALIZE THE ALPHA AND BETA ARRAYS 

IF (HOOIN,2I.EQ.Cl Nl=N/2 
IE IMODIN, 2I.E0. II Nl=IN-ll/2 
DC 10 J=!,N1 
DO 10 1=1,2 

N N =N— K K 

P T A 1 1 , J1=?1CNN » 

Ril'd , JI=22(NN 1 
KK =KK ♦ I 
IT CONTINUE 
N 2 =N 1 ♦ 1 
PT All ,N2> =1 . 

PT M2,N2» =1. 

IFINN.EO.il RTAC?,N2)=C. 

IF <NN.E0.2I RT All ,N21 =21 I 11 
IF INN. FQ. I I GO TO 20 
RTFC1 ,N21 =22 l 1 1 
PTE(2,N21 =0. 

20 CONTINUE 

CALCULATE THE ALPHA AND BETA ARRAYS 

DO 30 T=1 , N 
DO 30 J=1,N1 

AL FHA I I J=RTA ( I ,11 /RT A( 1*1 ,11 
BETAII J=RTeiI,ll/RTAII«l,ll 

RTMI«2,JI=RTA II , J*1 >- ALPHA 1 1 1 *R TA 1 1 ♦ 1 , J* 1 1 
RTEII»2,JI=RTBII,J*U-BETA(I1*RTAII*1,J*1» 

3 C CONTINUE 
RE TURN 
END 
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SUBROUTINE IMPULS 


C 
C 

C PURPOSE 

C TO COMPUTE THE ENERGY ASSOCIATED WITH EACH 

C ROUTH APP R CXI M AN T ASSUMING THE SYSTEM IS 

C DRIVEN WITH AN IMPULSE 

C 

SUBROUTINE DIMPLSI AL PH A E T A ,RM S ,N > 

IMFLICIT DOUBLE PRECISION «A-H,0-Z» 

DIMENSION ALPHA II I t BFT A f 1 > # R MS I 1 J 

S U M -D • 0 
DC ID T-l 9 N 

SUMBE TAII )**2 >12. * ALPHA f I I MSUM 
IF ISUM .LT .0.3 > RETURN 
R M S II ) -DSC R T CSLKI 
1C COMINUE 

return 

ENT 


r subroutine mpttrc 

c 

C PURPOSE 

C TO FORM THE STATE SPACE POUTH APPROXIMATE 

C OF THE ORDER N HA T FROM THE ORIGINAL RECIPROCAL 

C SYSTEM AND THE ROUTH TRANSFORMATION MATRIX 

C 

SUBROUTINE D R T IRC I AL F A f X I t E * C » T , T2 * TP t P t R HT , CH T f W A f WB , WC , 

I N f M t M M t N H AT } 

IMFLICIT DOUBLE PRECISION (A-H,0-?J 

DIMENSION ALF A CNI f PIN f MI , Cf MM,N 1 ,BHT INHAT # M ) ,CHT CMM t NH AT J 
] t V A I N 9 N I T WB I N ,N I , WC lMH f N UT lN f N » f TZ CN f N> ,TP IN ,N ) 9 DIM I ,XT CM 

FO^M THE ROUTH TRANSFORMATION MATRIX TP 

I J ID 

DO 1 J-UN 
IF II J. EC . 1 ) 60 TO 6 
IF IXII JJ.NE.C.n GO TO s 
DO II 1-1 ,N 
WB CI,J >-TZ< I ,J } 

II WA IlfJJrT IT,JJ$D(JJ 
GO TO 1 
E DO 12 T-l 9 K 

WB (I t J)-T2(J , J^IJ 
WLi Cl , J + l ) -T7I T ,J> 

VAII t J>:TU,JMDUMTIl f j4ll*DCj*l) 

1? WA II tJ + I! - - T I I vJ)4DfJ+19+T(I * J ♦ 1 ) * P I J ) 

I J = 1 
GC TO 1 
t I J -0 
1 CONTINUE 

CALL D TMUL T I W A f WB f TR f N ,N f NI 
I JN --1 
DC 2 J-I 9 N 
I JK--I JK 
ALFD-l JK*ALFA I 
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DO 2 1=1, N 

2 TR!1 # J)=TR(I,J1*ALFD 
FORM CHAT 

CALL DMULT!C»TR f WC,HM t N,N) 

DO 3 1=1, MM 
DO 3 J = 1 , NHA T 
2 CHTII , J»=WCt I ,J> 

F I ND THE INVERSE OF THF TR MATRIX 
I J -0 

CALL DHTMPYfWA ,UB,1.D0 ,N,N,N» 

DO 7 J = 1 , N 

1FIIJ.FQ.1J GO TO 72 
IF IX1C J1.EQ.C.01 GO TO 7 
DO 71 1=1 ,N 

V D <I,Jl=WAtI,J«l» 

71 Vh (1,J+1)=WA(I ,J) 

I J = 1 

GC TO 7 

72 I J =0 

7 CONTINUE 

CALL DTMULTIW8 ,TZ , WA , N ,N,N) 

I JK=-1 
DO 8 J = 1 « N 
I JK=— I JK 

ALFD=ALFA C J1*I JK 

00 8 1=1, N 

8 W A 1 1 , J ? -W ATI ,J 1*ALFD 

CALL DMUL T IWA,TZ,WP,N,N,N1 

I D C T = fc 

CALL OLNV1F l MB ,N,N ,WA ,TDGT ,UC ,IEP ) 
CALL DMULT!TZ,l<A,WB,N,N,N> 

FOFM BHAT 

CALL DHULTJWB ,e,WA ,N,N ,N1 
DO 4 I=1,NHAT 
DO 4 J = 1 , H 
4 BHT(I,JI=UA(I,J) 

RETURN 

ENC 


SUEROUTINE rthcnv 
PURPOSE 

TO COMPUTE THE COEFFICIENTS OF THr ROUTH 
CONVERGENT APPROXlKANTS 

SUEROUTINE DTHCNVI ALPHA.PET A ,NMMTX,DNMTX ,N,NMAX 1 
IMPLICIT DOUBLF PRECISION CA-H,0“Z» 

DOLBLE PRECISION NMMTX 

DIMENSION ALPHA! 1 I ,BET A ( 1 1 ,NMMT X I NM A X , 1 ),DNMTX!N*AX,1 J 
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DO 1C 1=1,NMAX 

ro 1C J-l ,NM ax 

KM 77 X ( 7 f J ) ~G • 

OK MX 1 1 f J)rr? . 

1C CC-M1NUE 

THE DENOMINATOR COP VE R GF MS ART 7 N PKMTX 
THE NUMERATOR CON VERGE N'T S ARE IN N M K T X 

PG ?C T=l f N 
PNHX4 1 ,7 i-l .a 
2C COMINUE 

rNMX<? t M=ALPEAU) 

DNMXf ?•? l-ALPHAC ? ) 
r N M X f 3 t r JrALF^AU 
N M M X C ? , 1 } - P E T A ( i > 

NM *T X <?*?>- r> ETAI?) 

n mmx< ? t? > - alpmi 2 **r r iA<n 

PC 30 I-3 9 * 

0 C 3C J~? t N 

CWMX4 J 9 I )“0NM7XCJ-lfI -1 )* ALPHA IT M PMM 7 X C J f I -2 ) 
NMMTXC J f J ) =MMM 7Xf J-l f I -1 I^ALPHA 11 1 *NMMTX< J» T ’2 1 
1MJ.EC.:) NMM 1XC J 9 1 ! = KMPTX < J 9 1 ) «bT T A C T I 
5 C COMINUE 
no Ir?,K 
DO ^0 J = 1 ,N 

PN MX I 7-1 f j) rOAMT X Cl , J 1 
N MMTX 4 7-1 *J)-NMMTX II »J ) 

ML CONTINUE 

ON MX <N,N*1 > =0 . 

NMMTX I N # N - 1 ) -G . 

P P CD - 1 . 

DO 5C 1~1 f N 
PRCD=PPOD^ALPHACI I 
5 C COMINUE 

DNMTXI N,N > -PROP 

NMMTX IN ,N )rpE 7 All > *PR 0 T / A L PH A ( I 1 
R T TURN 

e rr 


p SUBROUTINE DT Y PEI 

C 

c pi'pposr 

C TO FIND THf K 7 H ORDER TYPE I DOMINANT MOOT 

C APPROXIMATION OF AN NTH ORDER SYSTEM 

C 

C THF ORIGINAL SYSTEM 

C 2P0T r T^R^INVITI^Z + T*G*U 

C Y = H*INV<T)*? 

C 

SUBROUTINE D T Y FLI < A , B t C , T t T INV f R 1 1 , G »H t DK K 1 9 DKK2 9 DNN1 ,N,M f L,KI 
IMPLICIT DOUBLE PRECI S ION IA -H.O-Z ) 

DIMENSION AlK v M) v B(KtK>tCCL*KlvTCNfNl v PllCK t Kl«GIN 9 MI 9 H<LvNI 
1 «PKK 1 CK ,K ) 9 DKK2IH V K > f DNNl (N,N) T TINV4N f N) 

1 F IM.GT.K 1 WRI TEI fc ,100 Cl 
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irOC FORMAT (1HD ,1 DX ,'THt CONTPOL DIMENSION IS GREATER THAN THE* 
1 , IX . ’REDUCED ORDER DIMENSION IN DTVPEI’I 

I DCT-6 

FIND THE INVERSE OF THE K X K PLOCK OF T, Til 

DO 1 I = 1 , K 
DC 1 J = 1,K 
1 DKKMlfJl =T(I ,21 

CALL DMTMPVCDKKl,DKK? t l.DO.K ,K,K » 

CALL DLNVlF I 0 K K2 » K »K,A , 1 DGT , B , I E P I 

FIND THE CONTROL MATRIX P 


C 

c 

c 


CALL DMULT ITIN V,G , DNN 1 ,N , M ,N I 
DO 2 1=1, K 

DC 2 J=1,M 

D K t<2 I I , J) =ONN 1 II , J) 

CALL DMULTIDKK 1,DKK2,P ,K,M,K 1 

FIND THE OUTPUT M A TR I Y C 

CALL D M UL T IT , A ,CNN I ,N , K,K I 
CALL DMULTIH.DNN1 ,C,L,K,N» 

FIND THE SIMILARITY TRANSFORM T 1 1 *R 1 1 * I NV (Till 

CALL f; MULT IP 11 , A , D KK? , K ,K ,K ) 

CALL DMUL TIDKK 1,DKK2, A t K,K,K) 

R L TURN 
' U D 
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TABLE I. - COEFFICIENTS FOR TWO 16th-ORDER ENGINE 


TABLE H. - IMPULSE RESPONSE ENERGY RATIOS 


TRANSFER FUNCTIONS WITH FUEL FLOW INPUT 


Coefficient, 

Coefficient 

Compressor 

Turbine inlet 

i 

of characteristic 

speed, N c 

temperature, T t 


equation, 

(p = 1, q = 1) 

(p = 2, q = 1) 


a i 

Element of matrix of 



numerator coefficients B-, 


0 

1.0000 


1 

1. 0638xl0 3 

0.1140 

2 

3. 7805x10® 

4.2916X10 1 

3 

6.6915x10® 

-4.2414x10® 

4 

7.0216X10 10 

-2.2082x10® 

5 

4.7772X10 12 

-8.2163X10 7 

6 

2.2150X10 14 

2.2542X10 10 

7 

7.1949X10 15 

2.8702X10 12 

8 

1.6584X10 17 

1.5789X10 14 

9 

2.7151X10 18 

5.0230X10 15 

10 

3.1255X10 19 

1.0028X10 17 

11 

2.4739X10 20 

1.2850X10 1 ® 

12 

1.2973X10 21 

1.0446X10 19 

13 

4.2578X10 21 

5.1525X10 19 

14 

7.9983xl0 21 

1.4079X10 20 

15 

7.4304xl0 21 

1.7841X10 20 

16 

2.4119X10 21 

7. 4230X10 19 


5. 7270 

5. 8161x10® 

1. 8835x10® 

2.9209x10® 

2.6079X10 10 

1. 4714X10 12 

5. 5192X10 13 

1.4142X10 15 

2. 5027X10 16 

3.0507X10 17 

2. 5208X10 18 

1.3685X10 19 

IQ 

4.6358X10 13 
9. 0249X10 19 
8. 7170X10 19 
2. 9263xl0 19 


FOR INCREASING APPROXIMANT ORDER 


[Fuel flow input.] 


Routh 

approximant 

order, 

Compressor speed, N Q 
(p = 1, q = 1) 

Turbine inlet 
temperature, 
(p = 2, q = 1) 

k 

Impulse response energy ratio 

1 

0.41 

0.007 

2 

.72 

.03 

3 

.89 

.07 

4 

.94 

.14 

5 

.941 

.23 

6 

.95 

.36 

7 

.96 

.51 

8 

.98 

.67 

9 

.99 

.81 

10 

.997 

.91 

11 

.998 

.97 

12 

.999 

.99 

13 

.9992 

.998 

14 

.9998 

.999 

15 

.9999 

.9999 

16 

1.0000 

1.0000 



TABLE m. - ROUTH APPROXIMANT COEFFICIENTS FOR TWO 


TRANSFER FUNCTIONS 


Coefficient, 

i 

Compressor speed, 

N c (p = 1, q = 1) 

Turbine inlet temperature, T t 
(p = 2, q = 1) 

Routh approximant coefficients 


c. 

d. 

C- 

d i 


l 

i 

1 

1 

r\ 

i nnnn 


i nnnn 


u 

I , uuuu 




1 

2.6132 

0. 43978 

4. 7092 

3.6223 

2 

2.5702 

.61711 

1.0446X10 3 

9. 3469x10* 

3 

.8343 

.25676 

1.3858X10 4 

1.3400X10 3 

4 



1. 1722X10 5 

1.1906X10 4 

5 



6.3304X10 5 

6.6700X10 4 

6 



2. 1019X10 6 

2.2877xlO S 

7 



3. 9634X10 6 

4. 4718X10 5 

8 



3.6848X10 6 

4. 3229xlO S 

9 



1.1961X10 6 

1.4512x10 S 


TABLE IV. - COMPARISON OF POLE 
LOCATIONS FOR EXACT (16th ORDER) 
AND APPROXIMATE (3rd AND 9th ORDER) 


TRANSFER FUNCTIONS 


Poles of transfer functions 

16th Order 

3rd-Order 

9th-Order 


Routh 

Routh 


approximant 

approximant 

-0.648 

-0.635 

-0.648 

-1.91 

-0. 989±0. 580j 

-1.91 

-2.62 


-2.62 

-6. 71±1. 31 j 


-6. 50±1. 16j 

-17. 8±4. 80j 


-7.64±3.29j 

-18.2 

-21. 6±1. 55j 
-38.7 
-47. 1 
-50.6 
-59.2 
-175 
-577 


-6. 82±8. 73j 




TABLE V. - MINIMUM TRANSFER 
FUNCTION ORDER FOR 81-PERCENT 


IMPULSE RESPONSE ENERGY RATIO 


Transfer 

function 

output 

Transfer function input 

u i 

U 2 

U 3 

u 4 

u 

5 

x i 

3 

3 

4 

5 

HI 

x 2 

3 

3 

5 

4 


1 

x 3 

6 

12 

13 

12 

13 

x 4 

11 

12 

16 

11 

12 

x 5 

4 

5 

6 

4 

4 

x 6 

3 

6 

8 

3 



x 7 

3 

8 

7 

5 



x 8 

4 

5 

7 

7 



x 9 

8 

8 

7 

8 

8 

x 10 

4 

2 

2 

2 

A 


X 11 

9 

9 

7 

9 

9 

x 12 

8 

8 

7 

8 

8 

x 13 

7 

3 

5 

4 

6 

x 14 

7 

8 

7 

7 

6 


x 15 

4 

7 

7 

5 

4 


x 16 

6 

10 

9 

11 

10 
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TABLE VI. - BOOST- PUMP PRESSURE REGULATOR TRANSFER 


FUNCTION FOR PRESSURE ERROR MEASUREMENT INPUT P(s) 
TO ACTUATOR PISTON- COMMANDED VELOCITY OUTPUT Q(s) 

Q(s) _ -261. 31 s 42 - 9821. 8 s 41 - 1 39360 s 40 + . . . + 0. 34557x10® 

P(s) s 43 + 530. 64 s 42 + 1239. 2 s 41 + 17914 s 40 + . . . + 0. 33972X10 8 


Power of s 

Denominator 

coefficient 

Numerator 

coefficient 

42 

0. 5306 4xl0 2 

0. 26131X10 3 

41 

. 12392X10 4 

-.98218X10 4 

40 

. 17914x10® 

-. 13936x10® 

39 

.17766X10® 

-. 13898X10 7 

38 

.14079X10 7 

-.10524X10 8 

37 

. 91649X10 7 

-. 68170x10® 

36 

. 52194x10 s 

-. 37098X10 9 

35 

. 25635X10 9 

- . 18063X10 10 

34 

. 1146 7xl0 10 

- . 78084X10 10 

33 

. 45559X10 10 

-.30434X10 11 

32 

.16623X10 11 

-.10919X10 12 

31 

. 55692X10 11 

-. 35347X10 12 

30 

. 16873X10 12 

-. 10737X10 13 

29 

. 48960X10 12 

-.29612X10 13 

28 

.12415X10 13 

-.76725X10 13 

27 

. 31769X10 13 

-.18374X10 14 

26 

. 67590X10 13 

-.40591X10 14 

25 

. 15437X10 14 

-.85715X10 14 

24 

. 27519X10 14 

-.16043X10 15 

23 

. 56516X10 14 

-.30264X10 15 

22 

. 84066X10 14 

-.47440X10 15 

21 

.15582X10 15 

-.80866X10 15 

20 

.19203X10 15 

-.10439X10 16 

19 

.32126X10 15 

-.16249X10 16 

18 

.32473X10 15 

-.16879X10 16 

17 

.48853X10 15 

-.24235X10 16 

16 

.39945X10 15 

-.19629X10 16 

15 

.53594X10 15 

-.26269X10 16 

14 

.34801X10 15 

-.15896X10 16 

13 

.41055X10 15 

-.20056X10 16 

12 

.20670X10 15 

-.85553X10 15 

11 

.20970X10 15 

-.10314X10 16 

10 

. 79450X10 14 

-.28689X10 15 

9 

. 67087X10 14 

-.33582X10 15 

8 

. 18474X10 14 

-.55134X10 14 

7 

.12418X10 14 

-.63866X10 14 

6 

.23847X10 13 

-.54308X10 13 

5 

.11997X10 13 

-.63912X10 13 

4 

.15079X10 12 

-.22981X10 12 

3 

. 51938X10 11 

-.28645X10 12 

2 

. 39472X10 10 

-.18961X10 10 

1 

. 66059X10 9 

-.37172X10 10 

0 

. 33972X10 8 

. 34557X10 8 




TABLE Vn. - POLE COMPARISON BETWEEN APPROXIMANT AND ORIGINAL SYSTEMS 
(a) Approximant system (b) Original system 


Root 

Real part 
of root 

Imaginary 
part of root 

1 

-0. 57680x10"* 

0 

2 

-,17846xl0" 2 

.14393 

3 

17846 xlO' 2 

-.14393 

4 

-.29297X10' 2 

.29631 

5 

29297X10’ 2 

-.29631 

6 

-.51814x10"* 

. 34438 

7 

51814x10’* 

-.34438 

8 

-.57957X10" 2 

. 52760 

9 

-.57957X10" 2 

-.52760 

10 

43869x10* 

.69008 

11 

-.43869x10"* 

-.69008 

12 

-.37413x10"* 

-.10277x10* 

13 

-.37413x10'* 

-.10277X10* 

14 

-.35326x10"* 

.13448x10* 

15 

-.35326x10'* 

-.13448X10* 

16 

-.23819x10'* 

.14622x10* 

17 

-.23819x10'* 

-. 14622x10* 

18 

-.37308x10'* 

.17055x10* 

19 

-.37308x10"* 

-.17055x10* 

20 

-.38127x10"* 

.22899x10* 

21 

-.38127X10'* 

-.22899x10* 

22 

-.42813x10* 

0 

23 

-.55877x10* 

.43007x10* 

24 

-. 55877x10* 

-.43007x10* 

25 

-.44588x10* 

.60078x10* 

26 

-.44588x10* 

-.60078x10* 


Root 

Real part 
of root 

Imaginary 
part of root 

Root 

Real part 
of root 

Imaginary 
part of root 

1 

-0.57680x10"* 

0 

23 

-0.12750x10'* 

-0.17253x10* 

2 

-. 17846X10" 2 

. 14393 

24 

-.14813X10'* 

.19714x10* 

3 

-.17846X10" 2 

-.14393 

25 

-.14813x10'* 

-.19714x10* 

4 

-. 29297X10" 2 

.29631 

26 

-.11560X10"* 

. 22272x10* 

5 

-.29297X10" 2 

-.29631 

27 

-.11560x10'* 

-.22272x10* 

6 

-.51814x10"* 

. 34438 

28 

-.87273X10' 2 

. 24548X10 1 

7 

-.51814X10'* 

-.34438 

29 

-. 87273X10' 2 

-.24548x10* 

8 

-.57957x10" 2 

. 52760 

30 

-.64893X10' 2 

. 26466x10* 

9 

-. 57957X10" 2 

-.52760 


31 

-.64893X10' 2 

-.26466X10* 

10 

-.43869x10'* 

.69008 


32 

-.48184X10' 2 

. 27998X10* 

11 

-.43869x10'* 

-.69008 


33 

-.48184X10' 2 

-.27998x10* 

12 

-.37413x10'* 

. 10277x10* 


34 

-.36949x10' 2 

.29114x10* 

13 

-.37413X10'* 

-.10277X10* 


35 

-. 36949X10' 2 

-.29114x10* 

14 

-.32909X10"* 

. 13445X10* 


36 

-.30218X10' 2 

. 29785x10* 

15 

-.32909x10'* 

-.13445x10* 


37 

-.30218X10' 2 

-.29785x10* 

16 

-.14112x10'* 

. 14142x10* 


38 

-. 44449x10* 

. 44417x10* 

17 

-.14112x10'* 

-.14142x10* 


39 

-. 44449x10* 

-. 44417x10* 

18 

-.12075x10"* 

.15131x10* 


40 

-.91623x10* 

. 90855x10* 

19 

-.12075x10'* 

-.15131x10* 


41 

-.91623x10* 

-.90855x10* 

20 

-.18582x10"* 

. 16471x10* 


42 

-.42910x10* 

0 

21 

-. 18582x10'* 

-.16471x10* 


43 

-. 20926X10 2 

0 

22 

-.12750x10'* 

.17253x10* 






to 

CO 



TABLE Vin. - ZERO COMPARISON BETWEEN APPROXIMANT AND ORIGINAL SYSTEMS 


(a) Approximant system (b) Original system 


Root 

Real part 
of root 

Imaginary 
part of root 

1 

0. 91929x10' 2 

0 

2 

-.35351X10" 2 

. 14517 

3 

-. 35351X10" 2 

-.14517 

4 

-.29437X10' 2 

.29626 

5 

-.29437X10" 2 

-.29626 

6 

.92709X10' 2 

.35815 

7 

. 92709X10’ 2 

-.35815 

8 

-.72729x10' 2 

.52853 

9 

-.72729X10' 2 

-.52853 

10 

. 6 7278x10" 2 

.71196 

11 

. 67278x10 

-.71196 

12 

. 26794xl0~ 2 

. 10536X10 1 

13 

.26794x10' 2 

- . 10536x10* 

14 

-.57825X10" 2 

. 13715x10* 

15 

-.5782 5x10' 2 

-.13715X10 1 

16 

-.6 731 5x10' 1 

.14759X10 1 

17 

-.67315X10' 1 

-.14759X10 1 

18 

-. 22702x10' 1 

.17112X10 1 

19 

-.22702X10' 1 

-.17U2X10 1 

20 

.11499X10' 1 

.23645X10 1 

21 

. 11499X10" 1 

-.23645X10 1 

22 

-.32228X10 1 

.37158 

23 

-. 32228X10 1 

-.37158 

24 

-.30949X10 1 

-.40104X10 1 

25 

-. 30949X10 1 

.40104X10 1 


Root 

Real part 
of root 

Imaginary 
part of root 


Root 

Real part 
of root 

Imaginary 
part of root 

1 

0.91929X10' 2 

0 


22 

0.69556X10" 3 

0. 19874X10 1 

2 

-.29437X10' 2 

.29626 


23 

.69556X10' 3 

-. 19874X10 1 

3 

-. 29437X10' 2 

-.29626 


24 

-.15648X10' 2 

. 22398X10 1 

4 

.92709X10' 2 

.35815 


25 

-.15648X1 O' 2 

-. 22398X10 1 

5 

.92709X10" 2 

-.35815 


26 

-.26061X10" 2 

. 24639X10 1 

6 

-.72729X10' 2 

. 52853 


27 

-.26061X10' 2 

-.24639X10 1 

7 

-.72729X10' 2 

-.52853 


28 

-.29910X10' 2 

. 26526X10 1 

8 

. 67278X10” 2 

.71196 


29 

-. 29910x1 O' 2 

-.26526X10 1 

9 

.67278X10“ 2 

-.71196 


30 

-. 30101x1 0' 2 

. 28034x10* 

10 

.27081X10" 2 

. 10536X10 1 


31 

-. 30101x1 O' 2 

-.28034X10 1 

11 

. 27081X10" 2 

-. 10536 X10 1 


32 

-.29153X10" 2 

. 29131X10 1 

12 

-.66105X10" 2 

. 13790 X10 1 


33 

-.29153X10" 2 

-. 29131X10 1 

13 

-.66105x10' 2 

-. 13790 X10 1 


34 

-. 28165X10" 2 

. 29789X10 1 

14 

-.14819X10' 1 

.14161X10 1 


35 

-.28165X10" 2 

-. 29789X10 1 

15 

-.1 4819x10' 1 

-. 14161 X10 1 


36 

-. 35351X10' 2 

.14517 

16 

-. 51392x10' 1 

.15083X10 1 


37 

-.35351X10' 2 

-.14517 

17 

-.51392X10' 1 

-. 15083 xlO 1 


38 

-. 31477X10 1 

0 

18 

-.23310X10" 1 

. 16573 X10 1 


39 

- . 39426X10 1 

0 

19 

-.23310X10' 1 

-.16 573 X10 1 


40 

- . 46850X10 1 

. 43838X10 1 

20 

. 25798X10' 2 

.173 52 X10 1 


41 

-.46850X10 1 

-. 43838X10 1 

21 

. 25798X10' 2 

-.17352 X10 1 


42 

- . 20927X10 2 

0 





TABLE IX. - 26th-ORDER 


ROUTH APPROXIMANT 



Coefficient, 

Routh approximant coefficients 

i 

d i 

c i 

1 

-0.205702X10 3 

0. 249887X10 2 

2 

-.262997X10 4 

. 320663X10 3 

3 

-.190677X10 5 

. 248369X10 4 

4 

-.884145X10 5 

. 129942X10 5 

5 

-.315100X10 6 

. 485403xl0 S 

6 

-. 954589X10 6 

.151935X10 6 

7 

-.228406X10 7 

. 398134xlO S 

8 

-.481475X10 7 

.814672x10® 

9 

-. 887488xl0 7 

. 168262X10 7 

10 

-.129769X10 8 

. 234334xl0 7 

11 

-,199129x10 s 

. 398045X10 7 

12 

-.195556xlO S 

. 383160X10 7 

13 

-.263469X10 8 

. 540632X10 7 

14 

-.164054X10 8 

. 359295X10 7 

15 

-.202690X10 8 

. 416608X10 7 

16 

-. 739529X10 7 

.189273X10 7 

17 

- . 870168X10 7 

.175577X10 7 

18 

- . 170706X10 7 

. 539769X10 6 

19 

-.197216X10 7 

. 386031X10 6 

20 

-.186761X10 6 

. 789457x10 s 

21 

-.218839X10 6 

.411942x10 s 

22 

-.840811X10 4 

. 534361X10 4 

23 

-.103948X10 5 

. 188562X10 4 

24 

- . 726830X10 2 

. 144366X10 3 

25 

-.138053X10 3 

. 245338x10 s 

26 

. 128341X10 1 

. 126168X10 1 



TABLE X. - ENGINE 


[n =16, 1 = 2, 


(a) A 


-4.32800 

0.218859 

0.200384 

3.95571 

-2.98299 

-0.144833 

0. 769294X10" 1 

-0. 144819 

-.344744 

-5.64300 

3.72188 

-1.80121 

-1.40020 

1. 56028 

.120495 

-.299362 

27. 8479 

208.042 

-165.000 

-1.18466 

115. 858 

-165.738 

-10.7525 

20. 1663 

53.8484 

-14.0784 

496.865 

-578.300 

42.6308 

-70.2873 

-8.77498 

-36.0530 

2.05891 

-4. 84795 

. 507217 

3. 76362 

-10. 0500 

3, 55273 

-.161294 

-.755865 

11.1442 

-.212853 

-. 177448x10" 1 

-.798490X10" 1 

2.06863 

-19.7900 

-.182496 

-. 582995x10" 1 

8.97363 

-.213997 

-. 478545x10"* 

-.465845X10" 1 

4.37978 

-. 818499x10" 1 

-20.4700 

-. 768059x10' 1 

-.657472 

7. 50337 

4.26387 

. 701969x10 1 

-2.85571 

15.0693 

. 303936 

-19.9700 

-. 347415x10”* 

-. 826121x10”* 

.272257X10' 1 

.844885X10" 2 

-.137080 

. 207504 

. 168996X10" 1 

3.86258 

-.461982X10" 2 

-.110436X10" 1 

.360582X10' 2 

.112675X10" 2 

-.182567X10' 1 

.2 77078x10" 1 

.22 53 53x10" 2 

.515019 

-4.09217 

-29.0787 

3.31853 

. 417099 

-17. 7117 

26. 1649 

1.91267 

18.4196 

-.143460 

-3.04565 

-10.3264 

10.2900 

-.607676 

. 883145 

.71 1020x10" 1 

7. 55591 

-. 570369" 1 

-1.21806 

-4.13034 

4. 11573 

-.242984 

. 353474 

. 287463x10" 1 

3. 02213 

-.908447 

-1,77558 

-5.04877 

1.29700 

3.39206 

.691218 

. 588312x10" 1 

3.72598 

-.163698 

. 507438 

-. 587604x10" 1 

-.624872X10" 2 

.307427 

-.446147 

19.6607 

. 242382 

-11.9372 

16. 8839 

-1.02879 

8. 06630 

10. 1496 

-14. 7933 

-1.11110 

-.656673 


(b) B matrix 


-0. 469024x10" 1 

-0.124675 

. 895387x10" 1 

-.118071 

5.92813 

10.0878 

33.9885 

5.82984 

2. 48005 

-6.59813 

.373669 

. 263890 

. 233509 

.262489 

-. 538716 

-.264328 

. 156463 

-.841 755x10“ 2 

.208467X10" 1 

-.112409X10' 2 

19.9275 

1 -1.63330 

. 750732 

- . 352531x10 

.299580 

-. 141273X10“ 1 

.391725 

.2 58106X10" 1 

. 253898x10" 1 

. 332296x10" 1 

1.24292 

.868843 


^or p - 16, C is the 16th- order identity matrix. 
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MODEL DATA 


q =2 a ] 

matrix 

0.252614 
2.24174 
- . 636953x10 
-79.0958 
.128376 
-.659112x10 
. 629645x10 
.386312X10 
-49.9900 
-5.99940 
. 244677 
-12.7924 
-5.11614 
-6.28623 
-.375007x10 
.913777 

(c)C T 

1.00000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0. 427346X10“ 1 

-1.19810 

1.82058 

0.213133 

-0.106835 

-0. 997036X10 2 

0. 132963X10"* 

.239841 

0.940377 

-.855296X10" 2 

100398x10“ 1 

-.104123X10 -1 

-.353279x10“* 

-.104131X10“* 

.101911 

44.6143 

.114667 

. 101938 

.114658 

2.24267 

. 101904 

-12.9188 

170.896 

-4.53242 

-4.53157 

-4.43444 

-3.61533 

-4. 53105 

-.358717 

-2.92106 

-.424440 

-.420034 

1. 47850 

.922059x10“* 

-.422479 

-.608226X10" 1 

-. 253476x10“ 1 

-. 595834x10" 1 

-.6211 19x10" 1 

595750x10”* 

380239xl0“ 2 

-.621158X10“* 

-.793110X10" 1 

-.365148x10“* 

793274x10“ 1 

-.818338X10- 1 

793363x10“* 

.894036 <10“ 1 

-.818498x10“* 

.341217x10“* 

.167423X10" 1 

. 341306x10“ 1 

. 367042x10" 1 

.341328x10“* 

-.193187x10“* 

.341318x10“* 

.751111X10' 2 

46.0276 

. 375602X10" 2 

, 751138xl0“ 2 

. 375590x10 2 

0 

. 751331X10 -2 

-.665700 

6.13815 

. 450678x10“ 3 

.101370X10" 2 

.450645xl0“ 3 

0 

. 101 380x1 0~ 2 

.222139 

-47.6500 

.219747 

.219719 

.218157 

-.111896 

. 219763 

-1.41104 

59.6425 

-50.0100 

.112333X10" 1 

.112245x10“* 

187131x10“* 

. 124796X10“* 

564075 

23.8539 

-18.0000 

-1.99600 

. 449175x10 2 

-.718453X10" 2 

.494061x10“ 2 

-.688976 

29.3366 

-3.33640 

-.359398 

-19. 7700 

-.213957x10“* 

. 117686X10“ 1 

-.249931X10" 2 

. 712218x10“ 1 

249937xl0“ 2 

-.374785X10" 2 

-.250001X10" 2 

-20.0000 

-. 37482 4x10“ 2 

-. 291566x10“ 1 

-5.14614 

-.147456 

-.164604 

39.2746 

10.9155 

-50. 1600 


matrix 


0 


1 . 00000 
0 
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TABLE XI. - INPUT, OUTPUT, AND STATE VECTOR DEFINITIONS 


(Linear model: x = Ax + Bu; y = Cx, where x are state variables, y are output 


variables, and u are input variables. Since 

C = 1 (identity matrix), y = x.] 

State variables 

Control variables 

x 1 - Fan speed 

u^ - Main combustor fuel flow 

x 2 - Compressor speed 

u 2 - Nozzle jet area 

x 3 - Compressor discharge pressure 

u 3 - Inlet guide vane position 

x 4 - Interturbine volume pressure 

u^ - High variable stator position 

x & - Augmentor pressure 

u 5 - Compressor bleed 

Xg - Fan inside-diameter discharge temperature 


x 7 - Duct temperature 


Xg - Compressor discharge temperature 


Xg - Combustor exit temperature - fast component 


x 10 ~ Combustor exit temperature - slow component 

i 

- Combustor exit total temperature 


x 12 “ Fan ^rbine inlet temperature - fast component 


x^ 3 - Fan turbine inlet temperature - slow component 


x^ - Fan turbine exit temperature 


x^g - Duct exit temperature 


x^g - Augmentor exit temperature 



TABLE XII. - COMPARISON OF 

COMPUTATIONAL ACCURACY 

OF FREQUENCY- AND TIME- 

DOMAIN FORMULATIONS 

Alpha 

From Routh 
table 

0. 32460 
1.1230 
2.2886 
3.9662 
6.2954 
9. 3914 
13. 377 
18.425 
24. 806 
33.033 
44. 121 
60.072 
85. 236 
131.79 
254.39 
806.00 


parameters 

From 

transformation 

t r at r 

0 . 

32460 

1 . 

1231 

2. 

2886 

3. 

9662 

6. 

2954 

9. 

3914 

13. 

378 

18. 

425 

24. 

806 

33 

033 

44 

121 

60. 

.072 

85. 

.236 

131, 

.79 

254, 

. 39 

806 

.00 
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Figure 1. - Step response comparisons for normalized compressor 
speed and turbine inlet temperature outputs for fuel flow input - 
frequency-domain approach. 


Phase, deg Magnitude, dB 



1 2 4 6 8 10 20 40 60 100 

Frequency, Hz 

Figure 2. - Frequency response comparison. 
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